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(N; 

T-H , The densest packings of N unit squares in a torus are studied using analytical methods as well 

■ as simulated annealing. A rich array of dense packing solutions are found: density-one packings 
04 ■ when A'^ is the sum of two square integers; a family of "gapped bricklayer" Bravais lattice solutions 

\^ _ with density N/{N -\- 1); and some surprising non-Bravais lattice configurations, including lattices of 

^ . holes as well as a configuration for A = 23 in which not all squares share the same orientation. The 

■ entropy of some of these configurations and the frequency and orientation of density-one solutions 
' as A — >■ cxj are discussed. 

(^'. 

^ ; I. INTRODUCTION 

D. 

^ . Understanding the dense packings of hard particles has yielded essential insights into the structure of materials 
I ■ [IHll , granular media 0, Q , number theory 0] , biology d, Hj , and computer science [HI • This understanding 
, has been hard won: hundreds of years can elapse between a conjecture and its proof. This is best exemplified by 
sphere packing, for which a proof Kepler's conjecture had not been found until 1998 p^ . 

Recent experimental advances have allowed the development of (nearly) hard colloids that, for entropic reasons, 
"j^ ■ manifest the densest sphere packing [l^] . The attempt to obtain a deep understanding of liquid crystal mesophases 
' has prompted study of the densepacking of anisotropic particles. This has led to recent explorations of the packing 
7^ '. of ellipsoids ii3,[i|l, polyhedra [11, [13, and polygons [18]. 

■ One of the simplest regular polygons one can pack in two dimensions is the square. On the plane, the densest 
Ch , packing is, in this case, trivial - a square lattice of squares. Monte Carlo simulations of squares at finite pressure, 
^ ■ however, have also found a tetratic phase [1& . i2Q] . and experiments with hard colloidal squares have found, rather 

than the tetratic phase, a rhombic crystal having a different symmetry than the square ^l|. Even the dense packing 
of a finite number of squares can be more complicated than naive considerations would indicate - in fact, the problem 
] of packing squares in another square has been shown to be NP-hard [2^ . Indeed, the densest known packings can 
be quite complex [2^, [24| when the number of squares is not a perfect square integer. Higher packing densities than 
that of a simple, square lattice with vacancies can be achieved through configurations in which some of the squares 
^y!^ . are rotated and shifted with respect to the square lattice 
If^ ' What is the effect of an external potential on packings of hard objects? One example of an external potential is 
. fixed boundary conditions. Considerable effort has been devoted to understanding the densest packings of squares 

■ and circles in various domains. Another example with greater physical application is an external periodic potential. 
] If this potential is strong enough, the packing of the hard objects may be forced to be adopt the periodicity of 

^-H ■ the potential, and new packings are expected to arise. A limiting case of an imposed periodic potential is periodic 
^ ' boundary conditions. In this paper we explore the densest configurations of hard squares in a torus - that is, inside a 
. ^ . larger square with periodic boundary conditions. Even with the additional translation symmetry afforded by packing 

■ squares in a torus rather than in a square, the resulting dense packings in the torus can still be far from simple. Our 
. results may have experimental relevance for hard square colloidal particles in a periodic potential imposed either by 

. . . ■ a substrate or an optical lattice. 

As in many other mathematical packing problems, the strategy here is to search for the smallest area torus that can 
accommodate a fixed number of squares N . We use a combination of analytic and Monte Carlo simulated annealing 
techniques to accomplish this, and our results can be summarized as follows: we find that whenever N can be expressed 
as the sum of two square integers - N — nl -\- n2 - the densest possible configuration is a density-one packing with 
squares arranged in rows that are oriented at an angle of taii~^{n2/ni) relative to the underlying torus. For other 
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TV, we find a surprisingly rich collection of dense packing structures. For N = 6,11,14, and 27, we believe that the 
densest possible packing is a commensurate Bravais Lattice packing with density N/{N + 1) and resembles a bricklayer 
pattern with periodic gaps. For N — 12,21,22 and 23, we find that the densest configurations are non-Bravais lattice 
packings, including both regular lattices of holes and of rotated squares. These results are summarized the Table 
within Section II. 

In Section II, we present a summary analysis of the various structures we found for N up to 27, including both 
commensurate Bravais lattice solutions and non-Bravais lattice solutions. The packing motifs we found through 
analytic and numerical means are illustrated in this section via drawn figures as well as images generated by our 
numerical simulations. In Section III, we provide details of our numerical experiments for the hard square system. 
Section IV is a discussion of our results, including an analysis of the entropy of densest packings, and the rotational 
invariance of density-one packings as N goes to infinity. 



II. ANALYSIS OF PACKINGS AND NUMERICAL RESULTS 



In this section we give an analytic treatment of square packings on a torus. We first describe a set of solutions 
in which the squares lie on a Bravais lattice, and then turn to more complicated cases. The density-one solutions 
are optimal by construction, and all of the other solutions are conjectured to be optimal. The numerical results of 
Sec. mil guided us to the conjectured solutions, and the fact that long simulated annealing runs consistently produced 
these solutions gives us some confidence that they are optimal. Our conjectures for configurations containing N < 27 
squares are summarized in the Table. 



A. Commensurate Bravais Lattice Solutions 



Here we consider a class of Bravais lattice configurations that includes all of the density-one packings and other 
solutions we have found for N < 27. In all of these packings, the squares are lined up in rows; and for the purposes of 
this analysis, we assume that these rows are aligned along the x-axis. Thus, one of the primitive vectors of the lattice 
of squares is ai = x, and the second primitive vector is taken to have the form a2 = cx -|- dy, with —1 < c < 1 and 
\d\ > 1. Note that the primitive vectors of the torus, Ai and A2, need not be aligned with the primitive vectors of 
the squares. The requirement that the squares pack periodically on the torus is equivalent to saying that the lattice 
of squares is commensurate with the larger square lattice of the torus. That is, there exist integers ni, ri2, and 
such that the torus primitive vectors Ai and A2 are given by 



Ai 
A2 



riiai + n2a2 
ri^ai + 714 a2. 



In addition, we require that the torus primitive vectors are of equal length, 

|Ai| = |A2|, 

and orthogonal. 



These conditions are uniquely solved by 



Ai • A2 = 0. 



nin2 + n^Hi 



(1) 
(2) 
(3) 

(4) 



d = 



The number of squares N packed on the torus is the number of lattice points of the square lattice in a unit cell of the 
torus lattice 



N = \nini - n2n3|, 
and the areal density of the squares p is given by 

p = 7V/|Ai X A2| = l/|d|. 



(5) 
(6) 



3 



1. Density-one packings 

There are two classes of density-one packings. The first is the perfect square packing, for which c — 0, d = 1, 
ni — — VTV and n2 = = 0. This simple packing is shown in Fig. [T]for the case of = 9. Note that, on the 
torus, each of the ni rows (or columns, but not both) may be arbitrarily displaced relative to the other rows (columns) 
without disturbing the density of the packing or its periodicity; the perfect square packings thus have finite entropy. 



1 
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1 
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FIG. 1: An example of a perfect square packing with density one: A'^ = 9. 

There is a more general class of density-one packings, in which the lattice of squares may be tilted with respect 
to the primitive vectors of the torus. Setting d — 1 and c = in Eq. we find N — + n-l, rii = —714, and 
"■2 — ?^3j tfie square lattice being oriented at an angle of tan^^(n2/ni) relative to the torus lattice vectors. Clearly, 
these density-one, tilted square lattice solutions are optimal for all N that are sums of two square integers. Note that 
the perfect square solution corresponds to the special case n2 = ^3 = 0. Fig. [2] shows the case iV = 10 {rii — 3 and 
n2 = 1). 

We thus find that if N is the sum of two squares, there exists a density-one packing. The converse of this is also 
true: TV is a sum of two squares for all density-one packings of squares in the torus. To prove this, first note that 
every square in a density-one packing must have at least four other squares bordering it along a finite segment length, 
forcing all N squares to share the same orientation. Now consider three squares in mutual contact with each other - 
a configuration that must exist if the packing has no gaps. Two of those squares must be aligned in a row, as shown 
in Fig. [3l In order to eliminate gaps in the packing, these three squares define a set of rows that the entire packing 
must respect. Note that periodic boundary conditions allows us to draw the torus vectors so that they begin on the 
corner of a square and end on the corresponding corner of another square. Thus, 712 and 714 are both integers. A right 
triangle can be constructed with n2 as one side and Ai as its hypotenuse; another right triangle can be drawn with 
n4 as its base, and A2 as its hypotenuse. (See Fig. [3]). These triangles are identical by inspection. and 
are therfore each equal to + n1, via the Pythogorean theorem. 

2. Lattice packings with vacancies 

The simplest way to produce candidates for a densest packing for N = 122 + ni — k is to remove k squares from 
a density-one packing; indeed, our numerical results suggest that for several values of N, the densest packing is 
a density-one packing with one missing square. This is indicated in the Comment column in the Table using the 
notation — 1 or ri^ -|- — 1, depending on whether they are generated by removing 1 square from N a square 
integer, or a sum of two square integers respectively. As is demonstrated in Fig. 2] for the case of A^ = 2^ -|- 2^ — 1 = 7, 
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such vacancies allow for continuous displacement of other squares within a row, leading to a finite entropy for such 
configurations. Other examples include A'^ = 3 and 15. 

3. Bricklayer packings with gaps 

Next we consider Bravais lattice solutions that have density less than one - that is, packings with gap d — 1 > 0. 
Because these solutions have rows that are shifted relative to one another {c ^ 0), we call these "gapped bricklayer 
configurations." An example is shown in Fig. [SJ^a). Equations ^ allow us to ennumerate all gapped bricklayer 
configurations. Since the packing density of these configurations is {n\ + n\)/N , the highest packing density we can 
find within this class of configurations with density less than unity must have the form n\ + n\ = N — 1. This 
requires that N be one more than a sum of two squares. The first several are gapped bricklayer configurations are 
seen for N — 2,3,6,11,14,18,26, and 27. Based on the numerics we believe that for = 6, 11, 14 and 27, the gapped 
bricklayer packing is the densest possible packing. These are indicated in the Table with the abbreviation "GB" in 
the Comment column. Associated (non-unique) lattice vectors are shown in the rightmost columns of the Table for 
the GB packings. There are also gapped bricklayer solutions for density (N — 2)/N when N is two more than a sum 
of two squares, though we have not found any candidate densest packing solutions of this form for N < 27. Unlike 
the bravais lattice packings with density one, different rows of the gapped bricklayer solutions have a fixed shift given 
by c = —{nin2 + njn4)/ln\ + nl), where the denominator is the closest sum of two squares below N . 

B. Non-Bravais Lattice Packings 

Here, we consider special cases suggested by our numerical simulations that do not correspond to Bravais lattice 
packings. Note: if the Comment column of the Table simply repeats a value of N ^ it indicates a special case for which 
the squares are not on a Bravais lattice and for which there is not an obvious pattern that can be extrapolated easily 
to optimal packings for higher values of TV. 

1. Gapped bricklayer with domino bricks, N = 22 

The conjectured best packing for iV = 22 is shown in Fig. [6l This packing is in fact a gapped bricklayer configuration, 
except that the unit cell or brick is composed of two squares stacked in the y direction (the direction perpendicular 
to the rows). The configuration is otherwise identical to the = 11 gapped bricklayer and has density p = 10/11. 

2. Lattice o/ i x i holes, N = 12 and 23 

The conjectured best configurations for = 12 and 23 are shown in Figs. [7] and [3 respectively. In both cases the 
motif can be described as a lattice of ^ x ^ holes. The torus lattice vectors, Eq.[T]with ai = x and a2 — y, for A^ = 12 
are described by ni = n2 = —n^ = 714 = 5/2 and the torus lattice vectors for A^ = 23 are described by ni = = 9/2, 
and —712 = '1-3 = 2. It is straightforward to verify that these motifs are in fact packings on the torus and have the 
density N/ {n\ + n|) = N / {N + fc/4) where k is the number of holes in the unit cell. For A^ = 12, evidently k — 2 and 
for A^ = 23, fc = 5. 

3. Lattice of skew squares embedded m a square lattice, N = 21 

The conjectured densest packing for A^ = 21, shown in Fig. [HI does not follow any of the motifs described heretofore. 
The unit cell consists of a 4 x 4 square with motif of 5 squares attached to its side. This 5-square pattern is also the 
best packing of 5 squares in a square [13]. A simple calculation yields the density, p = 21/(4^ + (2 + l/\/2)^). This 
packing has one square per unit cell tilted at 45° relative to all other squares. This is the only example that we found 
for which not all of the squares in the motif are oriented in the same way. It was also the most difficult configuration 
for our simulated annealing algorithm to find. Figure [9] shows a typical simulation result, which clearly has not yet 
fully converged. 
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TABLE I: Exact and conjectured densest packed configurations of squares on a torus. Refer to tiie text for the meaning of tiie 
columns. 



N 


P 


Comment 


ni 


n2 




n4 


1 


1 


1^ 


1 








2 


1 


1^ + 1^ 


1 


1 






3 


3/4 = 0.75 


2^-1 


2 








4 


1 


2^ 


2 








5 


1 


2^ + 1^ 


2 


1 






6 


5/6 = 0.83 


GB 


2 


-1 


2 


2 


7 


7/8 = 0.875 


2^ + 2^-1 


2 


2 






8 


1 


2^ + 2^ 


2 


2 






9 


1 


3^ 


3 








10 


1 


3^ + 1^ 


3 


1 






11 


10/11 = 0.90 


GB 


3 


1 


-2 


3 


12 


24/25 = 0.96 


12 










13 


1 


3^+2^ 


3 


2 






14 


13/14 = 0.9285714 


GB 


4 


-2 


1 


3 


15 


15/16 = 0.9375 


4^-1 


4 








16 


1 




4 








17 


1 


42 + l2 


4 


1 






18 


1 


3^+3^ 


3 


3 






19 


19/20 = 0.95 


4^ + 2^-1 


4 


2 






20 


1 


4^+2^ 


4 


2 






21 


21/(4^ + (2 + 1/^2)^) = 0.900189 . . . 


21 










22 


10/11 = 0.90 


22 


3 


1 


-2 


3 


23 


92/97 = 0.94845 . . . 


23 










24 


24/25 = 0.96 


52-1 


5 








25 


1 


52 


5 








26 


1 


5^ + 1^ 


5 


1 






27 


26/27 = 0.962 


GB 


5 


1 


-2 


5 



C. Table of Results 

To summarize: the perfect square, sum of two squares and gapped bricklayer configurations cover most of the case 
we have found for N < 27. The Table gives densest packing configurations (if p = 1) and conjectured densest packing 
configurations (if p < 1), for each value of less than 28. The column "p" is the density of the configuration. The 
Comment column describes the type of lattice. For example, 3^ indicates a perfect square and 5^ — 1 indicates a perfect 
square with one square missing. Similarly 3^ + 1^ refers to the sum of two squares and "GB" stands for "gapped 
bricklayer." If a single number is in the Comment column it refers to one of the special cases discussed above. The 
four columns "ni, n2, n^, 714" are shown if the configuration of squares is itself a Bravais lattice; these integers are 
the coefficients of the lattice vectors of the torus, in terms of the lattice vectors of the squares as defined in Eq. ((T|). 
Only those columns needed to specify the lattice are filled in. 



III. NUMERICAL METHODS 

For all N < 27 squares on the flat torus, we searched for densest packings of N squares on the torus via Monte 
Carlo simulations in the NPT ensemble. Our approach was to employ a simulated annealing (SA) algorithm in which 
the system was taken from an initial, low-pressure, easy-to-equilibrate state to a final, high-pressure state, via an 
annealing schedule consisting of a series of steps in inverse pressure. Between each pressure increase a Metropolis 
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algorithm appropriate to the hard square NPT ensemble was used to equilibrate the system. Although SA quickly 
falls out of equilibrium at higher pressures as the energy landscape becomes rough, it appears to be an an effective 
algorithm for finding ground states of the system. 

The equilibration procedure we used in our simulated annealing algorithm was a Metropolis procedure consisting of 
three types of Monte Carlo moves: translations and rotations of individual squares, and changes in the volume of the 
entire system. At each step of the equilibration procedure, a square is selected at random; then, one of the three types 
of moves is selected at random, with probabilities .495, .495, and .01, for translation, rotation, and volume change, 
respectively. Once a square and a move type is selected, the move is attempted. If the move results in any overlaps 
among the N squares, the move is rejected. If it does not, then for translations or rotations, the move is accepted; 
for volume change dV, the move is accepted with probability Pacc — min[l, exp(— /JPdy)]. In practice, rather than 
changing the volume of the entire system, the periodic box in the simulations was kept at a constant size, and the sizes 
of the individual squares were all rescaled, in order to achieve the desired new volume. The equilibration procedure 
consists of s such Monte Carlo steps; in our simluations s was typically between 200 and 400 steps. 

We repeated this simulation 1,000 times and reported the highest density found among these runs. In order 
to determine the highest-density packing for N that did not correspond to a perfect square or a sum of squares, 
more extensive runs were conducted - in some cases, as long as 72 hours on a 2GIIz processor. The fact that 
significantly different initial configurations as well as different initial random seeds generated the same final, high- 
density configurations signalled that a good candidate for a densest packing of the system had been found. For all 
simulations, the pressure was initially set to l3P — .01, and was increased via constant steps in inverse pressure until 
a maximum pressure of pPmax = 3000 was reached. (During subsequent explorations of the phase behavior of the 
system, this pressure was later deemed excessive, but nevertheless produced reasonable results for the purposes of 
determing the ground state of the system.) All simulations were begun at an initial areal density p of 0.1, with a 
square array of unit squares. Before each equilibration procedure began, a trial run was conducted in which the 
maximum value of translations, rotations, and volume changes were independently optimized in order to achieve an 
acceptance ratio for each of 0.4. 

The majority of the computational work during the equilibration procedure consisted in checking for square overlaps. 
For this, we relied on a fast algorithm for detecting polygon overlaps by Alan Murta |.25i]. and an associated Python 
wrapper by Joerg Radler [2^ . Configurations were visualized using the VPython library [23| • 

IV. DISCUSSION 

In this paper, we have presented an analysis of the densest packing solutions for N unit squares in the torus. For 
N < 27, the majority of these packings are Bravais lattice solutions, manifested either as the "sum of two square 
integers", or "gapped bricklayer configurations" described above. A few were non-Bravais lattice solutions, such as 
those = 21 and iV = 23. In this section, we discuss the frequency and entropy of the various types of packings we 
found, and pose some questions for further study. 

We showed in Section fll A II that density-one packings are only possible for those values of N that are expressible 
as the sum of two squares. Though it appears that density-one packings are relatively common from the Table, in 
fact it is known that the frequency of numbers that are equal to the sum of two square integers scales as 1/VhiN for 
large N [28, J9j- Thus the frequency of density-one packings also vanishes with increasing N. 

Despite the relative scarcity of density-one packings, we argue that the packing density approaches one as A^ ^ oo. 
Suppose M is a sum of two square integers. Construct a new packing for N = M — k by removing k squares. The 
resulting packing density gives a lower bound of p = N/{N + k). Given a number of squares, N, however, determining 
k requires knowledge of the nearest density-one packing la rger t han N. For sufficiently large N, we can estimate that 
the distance to the next density-one packing is of order vhiiV larger than N. Therefore, an estimated lower bound 
on the packing density of N squares is 1 — Vln N /N for large N, which yields the result that the density approaches 
one asymptotically. This argument does not take into account fluctuations in the spacing of sums of two squares, and 
it would be interesting to find a mathematically rigorous asymptotic lower bound on the packing density. 

The main contributions to the entropy of the various Bravais and non-Bravais lattice configurations for N < 27 are 
readily assessed. For a Bravais lattice packing, there are no non-trivial continuous symmetries and thus there is no 
entropy. However, one can construct density-one packings from the Bravais lattice packings by shifting rows relative 
to one another. For example, if A is a perfect square then each row can be arbitrarily shifted and, ignoring overall 
shifts of the lattice, the entropy is proportional {Vn — 1). Note that either rows or columns may be shifted but not 
both for density-one packings. More generally, rows are free to shift unless the contraints of periodicity forbid it. 
Periodicity forbids two rows from shifting relative to one another if some linear combination of torus lattice vectors 
connects the rows. Suppose that rows are aligned in the x-direction. Then, two rows are constrained to their relative 
positions in the Bravais lattice configuration if there are integers a and b such that aAi -|- 6A2 connects the two rows. 
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Thus, rows separated by an2 + bn4 are locked (see Eq. [T]and Fig.[3](b)). Bezout's Lemma [so'] states that this integer 
hnear combination can be made equal to, but not smaller than, g, where g is the greatest common divisor of 712 and 
n4 (assuming that both 712 and 714 are nonzero). Thus, the set of rows can be divided into g groups, each of which can 
be arbitrarily shifted, and the entropy is proportional to g—1. Note that if I712I and |n4| are mutually prime, all rows 
are locked and this contribution to the entropy of the configuration vanishes. There are two other sources of entropy 
for packing related to Bravais lattice packings. The gapped bricklayer configurations (Section III A3|) allow squares 
within a row to shift perpendicular to the row axis (see Figures [S] and ^ , and this "poor workmanship" contribution 
to the entropy will be roughly proportional to the free volume of the configuration. The density-one configurations 
with vacancies (n2 + n1 — k), discussed in Section lll A 21 also possess finite entropy, since the hole(s) created by the 
k missing squares may be moved throughout the lattice, or split along a row; and for fc > 1, holes can appear in 
different rows. In contrast, the unusual, non-Bravais lattice packings of the = 12 and = 23 exhibit no entropy. 

The above entropies have implications for which configurations are most likely to be seen at finite pressure. For 
example, A^ = 25 admits of two classes of configuration: a density-one packing with rows aligned along the torus 
(A = 5^) or with rows oriented at an angle with respect to the torus (A^ = 3^ + 4^). As we have seen, the N — 5'^ 
has four rows that are free to slide while TV = 3^ -I- 4^ is a locked configuration with no entropy; this implies that the 
N = 5'^ packings will be much more likely to appear at finite pressure. 

Unlike squares packed into a square boundary, squares packed on a torus maintain rotational invariance in the 
thermodynamic limit. This can be seen as follows: one can see in Figure [5] that any N — nl + n^ packing (which are, 
as seen above, the only possible packings with density one) will orient the square lattice at an angle of tan~^ (77,2,1/711,1) 
relative to the torus lattice vectors. To take the thermodynamic limit with a particular square lattice orientation Q 
with respect to the underlying torus, it is sufficient to choose a particular subsequence of integers Ni = n\ ^ + n\ ^ such 
that tan"-'^ (712/711) — >■ Q. The thermodynamic limit of density-one packings on the torus thus preserves rotational 
symmetry. 

Our study of the densest packings of A^ unit squares in a torus has yielded definitive results for cases in which A^ is 
the sum of two square integers or is a perfect square, and strong conjectures for other values of A^ < 27. This work 
raises many interesting questions: How common are densest packings that have squares with different orientations, 
such as occurs for A^ = 21? Which motifs, if any, dominate for large A? Are the 1/2 x 1/2 motifs of A^ = 12 and 
A = 23 exhibited for other A^? What role do these dense packing configurations play in the thermodynamic phases 
observed in finite-temperature simulations and experiments with square colloids? 
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FIG. 2; An example of a packing for which TV is equal to the sum of two squares: A'^ = 10; all such packings are density one. 




(a) (b) 



FIG. 3: (a) Any three squares in mutual contact without gaps force two of the squares to define a row. (b) Diagram illustrating 
that A must be equal to the sum of two squares for all density-one packings of squares in the torus (see discussion in section 
UTAlli . 
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FIG. 5: Schematic (a) of a "gapped bricklayer" configuration, with density p = {N — 1)/N. Results of simulated annealing 
for A*' = 11 arc shown in (b) (rii = 3, n2 — l,n3 = — 2,n4 = 3). The finite entropy of this configuration is revealed by the 
displacements of the squares perpendicular to the rows. 




FIG. 6: Simulation results for N = 22: the conjectured best packing is a "gapped bricklayer with domino bricks" . 



FIG. 7: Simulation results for N = 12: a lattice of i x i holes. 




FIG. 8: Simulation results for N = 23: a lattice of | x i holes. 
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FIG. 9: Simulation results for — 21: a lattice of skew squares embedded in a square lattice. The 5-square pattern (which 
includes a skew square in its center) is the proved best packing of 5 squares in a square [13]. Note that the simulation results 
have not yet converged to the conjectured best packing. 



